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Abstract 

The  methodology  of  vector  and  matrix  algebra  has  been  used  to  simplify  the 
triangulation  procedures  required  to  determine  the  position,  motion,  and  growth  of 
luminescent  gas  clouds  injected  into  the  upper  atmosphere.  Generated  by  chemicals 
discharged  from  rockets  programed  to  effect  point  or  continuous  releases  at 
preselected  times,  such  clouds  have  proved  to  be  a  powerful  experimental  tool  for 
gathering  data  relevant  to  upper  atmosphere  phenomena  such  as  winds,  wind  shears, 
and  turbulent  transport  mechanisms,  and  the  production,  maintenance,  and  decay 
of  ionization  in  support  of  specific  Air  Force  requirements  in  satellite  operations, 
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Application  of  Vector  and  Matrix  Methods 
to  Triangulation  of  Chemical  Releases 
in  the  Upper  Atmosphere 


1.  INTRODUCTION 

Photographic  data  taken  from  ground  stations  for  the  purpose  of  determining 
the  position,  motion,  and  growth  of  luminescent  clouds  and  trails  have  for  a  number 
of  years  been  reduced  and  analyzed  by  various  procedures  (Refs.  2,5,6, 7,8, 12). 

The  artificial  clouds  and  trails  are  generated  by  discharging  chemicals  into  the 
upper  atmosphere  from  rockets  programed  to  effect  point  or  continuous  releases 
at  preselected  times.  They  have  proved  to  be  a  powerful  experimental  tool  for  the 
study  of  various  upper  atmosphere  phenomena,  including  winds,  wind  shears,  and 
turbulent  transport  mechanisms. 

The  triangulation  techniques  applied  to  the  data  differ  according  to  the  type  of 
release  whose  evolution  is  recorded  on  the  film  being  processed.  As  might  be 
anticipated,  the  analysis  becomes  simpler  for  discrete  puffs  possessing  structural 
features  that  can  be  clearly  identified  from  two  or  more  ground  stations.  A  more 
complex  mathematical  problem  is  presented  by  diffuse  clouds  and  continuous  trails 
with  few  or  no  distinctive  features  that  can  be  unambiguously  matched  on  the  film 
records  corresponding  to  views  of  the  release  from  different  triangulation  sites. 

The  present  report  describes  an  attempt  to  simplify  reduction  and  interpreta¬ 
tion  of  the  photographic  data  and  increase  accuracy  and  reliability  of  the  results 
by  (a)  devising  analytic  procedures  that  consistently  discard  superfluous  coordinate 
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transformations,  (b)  try  to  develop  formulas  that  are  symmetric  with  respect  to 
the  observation  sites,  (c)  establish  direct  analytic  criteria  for  determining  coor¬ 
dinates  of  the  cloud  features  under  observation,  and  (d)  use  redundancy  to  check 
the  consistency  of  previously  determined  results.  The  reduction  in  the  number  of 
coordinate  transformations  leads  to  a  conceptually  simpler  computational  scheme, 
and  as  direct  consequence  of  fewer  arithmetical  operations  having  to  be  performed 
on  the  original  data,  the  accuracy  of  the  results  is  degraded  to  a  lesser  degree. 

Mathematical  formulation  of  the  triangulation  problem  becomes  very  straight¬ 
forward  by  systematic  use  of  vector  algebra  and  analysis.  Recognition  of  the 
vector  character  of  the  positional  data  allows  many  of  the  expressions  to  be  written 
in  invariant  fashion,  and  the  coordinate  transformations  become  simple  rotations 
and  translations.  The  abstract  properties  of  the  corresponding  operators  are 
realized  by  well-known  matrices,  and  the  entire  scheme  can  be  cast  into  a  form 
that  can  be  generalized  with  minimum  effort.  The  resultant  expressions  are  ex¬ 
tremely  compact  and  can  often  be  readily  grouped  into  symmetric  combinations  in 
which  data  from  all  sites  play  equal  roles.  Aside  from  the  fact  that  there  is 
generally  no  a  priori  reason  for  assigning  a  preferential  position  to  a  given  site, 
an  advantage  of  such  a  scheme  is  that  special  or  singular  cases  and  the  concomitant 
lengthening  and  increased  complexity  of  the  computer  codes  needed  to  implement 
the  mathematical  formalism  are  avoided.  As  a  further  advantage,  well-known 
vector  identities  can  be  used  to  check  the  consistency  of  the  results. 

It  is  hoped  that  the  organization  of  the  material  will  allow  the  reader  to  follow 
development  of  the  thesis  without  undue  effort.  At  the  outset  we  recognize  that  die 
triangulation  problem  can  be  divided  into  independent  units.  Some  are  repeatedly 
used  during  analysis  of  a  batch  of  data.  Others  are  needed  only  to  fix  the  orienta¬ 
tion  of  the  cameras  and  are  not  recalled  unless  the  cameras  are  moved  during  the 
time  it  takes  to  photograph  the  entire  sequence  of  events  associated  with  the  deposi¬ 
tion,  growth,  and  decay  of  a  chemical  release.  Orientation  and  some  optical 
parameters  of  the  camera  can  be  accurately  determined  by  photographing  the  stellar 
background  shortly  before  the  rocket  and  its  chemical  payload  are  launched.  Stars 
can  be  readily  identified  and  their  positions  on  the  photographic  plate  can  be 
measured  very  precisely.  The  plane  of  the  photographic  plate  and  the  line  of  sight 
along  the  camera  optic  axis  serve  to  define  one  of  the  coordinate  systems  of  direct 
interest  to  us.*  The  other  systems  are  the  equatorial  system,  traditionally  used 


*In  the  discussion  presented  here  it  is  assumed  that  a  perfect  camera  is  used. 
Performance  of  actual  cameras  must  be  corrected  for  such  defects  as  shift  of 
center  of  frame  from  trace  of  optic  axis,  departure  of  film  plane  from  perpen¬ 
dicularity  to  optic  axis,  lens  distortion,  and  film  shrinkage.  The  analysis  remains 
unaltered  if  it  is  understood  that  appropriate  corrections  have  been  applied  to  all 
film  measurements.  For  procedures  to  evaluate  the  corrections,  see  Justus 
(1963). 
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by  astronomers  as  a  frame  of  reference  for  star  locations,  and  a  cartesian 
geocentric  system  with  its  origin  at  the  center  of  the  earth,  its  +z  axis  along  the 
lire  joining  the  origin  and  the  North  Pole,  its  +x  axis  coincident  with  the  line 
joining  the  origin  and  the  point  on  the  equator  intersected  by  the  prime  meridian 
(Greenwich),  and  its  +y  axis  completing  a  righthanded  orthogonal  triad.  In  view 
of  the  fact  that  our  observations  are  made  from  stations  riding  on  the  earth,  it  is 
advantageous  to  think  of  this  geocentric  system  as  a  preferred  frame  cf  reference 
•  and  to  regard  geocentric  coordinates  as  canonical  coordinates  for  all  events.  In 

particular,  this  will  simplify  the  conversion  from  equatorial  coordinates  as  well 
as  the  conversion  to  two  auxiliary  systems  we  shall  use  on  occasion:  an  azimuth- 
elevation— slant  range  system  and  the  geodetic  system  of  longitude,  latitude,  and 
altitude  above  surface  of  a  reference  ellipsoid. 

Section  2  deals  with  the  basic  coordinate  systems.  The  coordinate  trans¬ 
formations  are  developed  in  Sec.  3;  in  Sec.  4  they  are  applied  to  the  previously 
defined  systems.  At  this  point  we  have  sufficient  information  to  deal  with  the 
problem  of  determining  the  geocentric  coordinates  of  a  point  release  observed 
from  two  or  more  stations.  This  constitutes  the  bulk  of  Sec.  5. 

Section  6  describes  the  more  elaborate  procedure  that  is  needed  to  triangulate 
on  thin,  continuous  trails  of  lengths  measured  in  tens  of  kilometers.  The  most 
vexing  problem  associated  with  these  trails  is  the  difficulty  of  matching  correspond¬ 
ing  points  on  photographs  taken  at  different  tiroes  or  from  different  stations.  The 
procedure  we  follow  is  to  select  a  point  on  one  view  and  search  for  its  image  on  a 
second  view  after  the  trail  has  been  replaced  by  a  two-dimensional  curve  on  the 
plate  plane,  which  is  itself  represented  analytically  by  a  parametric  cubic  spline. 

The  introduction  of  parametric  representations  is  mandatory,  to  deal  with  multi¬ 
valuedness  created  by  selfintersections  or  pronounced  twisting  of  the  trail  and  its 
replacement  curve.  Since  very  little  can  be  found  in  the  literature  on  parametric 
spline* ,  the  general  procedure  needed  to  determine  the  sets  of  polynomial 
coefficients  is  described  in  Appendix  A. 

Section  7  deals  with  coordinate  transformations  between  geodetic  and  geocentric, 
geocentric  and  horizon,  and  camera  and  horizon,  systems. 

Since  all  geodetic  reference  ellipsoids  have  very  small  eccentricities,  it  was 
found  convenient  to  make  use  of  the  Lagrange  expansion.  A  statement  of  the  ex¬ 
pansion  in  its  full  generality  is  given  separately  in  Appendix  B.  Appendix  C  lists 
a  short  table  of  useful  constants. 

[ 

!  I 

2.  THE  COORDINATE  SYSTEMS 

Systems  of  reference  for  astronomical  purposes  are  constructed  by  selecting 
l  a  great  circle  on  the  celestial  sphere  and  one  point  on  that. circle.  Spherical 

l  ' 
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coordinates  are  then  introduced  to  specify  the  position  of  points  on  the  sphere: 
one  coordinate  measured  perpendicular  to  the  selected  great  circle  along  an 
auxiliary  great  circle;  the  second  measured  from  the  selected  point  to  the  point  of 
intersection  of  the  auxiliary  circle. 

The  fundamental  astronomical  reference  systems  are  based  on  the  celestial 
equator  and  the  ecliptic.  Angular  coordinates  in  these  planes  are  measured  from 
the  ascending  node  of  the  ecliptic  on  the  celestial  equator,  customarily  referred  to 
as  the  vernal  equinox,  the  first  point  of  Aries,  or  simply  the  equinox,  and  denoted 
by  T .  Only  the  first  of  these  systems  is  relevant  to  our  discussion. 

A  righthanded  cartesian  coordinate  system  with  its  origin  at  the  center  of  the 
earth  can  be  constructed  by  directing  the  +x  axis  toward  the  equinox,  the  +y  axis 
to  a  point  »/2  radians  to  the  east  in  the  plane  of  the  celestial  equator,  the  +  z  axis 
pointing  north  to  the  pole.  This  is  the  equatorial  system,  shewn  schematically 
in  Figure  1.  To  specify  the  location  of  very  distant  points  it  is  sufficient  to  give 
the  angular  coordinates  that  identify  a  direction— customarily  the  angles  a  and  6, 
which  denote  the  right  ascension  and  declination,  respectively.  Clearly,  a  and 
(»/ 2)-6  are  the  familiar  spherical  coordinates  corresponding  to  that  direction  in 


Figure  1.  Equatorial  Coordinate  System 


space.  If  necessary,  the  same  direction  can  be  identified  by  the  unit  vector  e 
having  the  cartesian  representation 
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e  «  x  (cos  a  cos  6)  +  y  (sin a  cos  6)  +  z  sin  6  . 

Any  point  along  the  line  OP  can  tnen  be  represented  by  the  vector 

A 

r  =  r  e. 


(1) 


(2) 


where  r  denotes  the  distance  from  O. 

The  declination  is  measured  in  degrees,  minutes,  and  seconds  along  the  great 
circle  through  the  observed  body  and  the  pole.  Right  ascension  is  measured  in 
hours,  minutes,  and  seconds  from  T  along  the  equator  in  the  west-east  direction, 
that  is,  opposite  to  the  apparent  rotation  of  the  celestial  sphere.  The  conversion 
to  the  equivalent  angular  measure  offers  to  difficulty. 

Since  observations  are  made  from  the  surface  of  the  earth,  topocentric 
coordinate  systems  with  origin  at  the  obser  -er  are  of  importance.  The  reduction 
from  earth-centered  to  topocentric  coordinates  depends  in  part  on  the  figure  of 
the  earth  and  will  be  considered  in  later  sections. 

For  our  purposes  it  'S  advantageous  to  inti  oduce  as  the  primary  geocentric 
system  a  righthanded  cartesian  system  fixed  in  the  earth,  with  the  +z  axis  pointing 
north,  and  the  equator  lying  in  the  xy  plane.  The  +x  axis  is  directed  toward  the 
intercept  of  the  equator  and  a  prime  meridian  (Greenwich,  for  example).  The 
+y  axis  completes  an  orthogonal  triad  (see  Figure  2..  The  context  will  be  sufficient 
to  distinguish  this  system  from  others  that  share  the  same  origin.  The  angular 
coordinates  of  this  geocentric  system  will  be  described  by  using  X  and  <t>’  to  denote 
the  longitude  and  geocentric  latitude,  respectively.  A  unit  vector  g,  which  in  the 
geocentric  system  can  be  written  in  component  form  as 


g  =  x  ( cos X  cos <t>')  +  y  (sin*-  cos  $')  +  z  sintf' 


(3) 


corresponds  to  the  direction  specified  by  the  coordinates  (X,  4>').  The  position 
vector 


r  =  rg 


(4) 


corresponds  to  a  point  along  this  direction,  r  units  away  from  the  origin. 

The  basic  topocentric  system  of  interest  to  us  is  determined  by  the  cameras 
that  record  the  evolution  of  the  chemical  releases. 
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Figure  2.  Geocentric  Coordinate  System 


The  +z  axis  of  the  coordinate  system  associated  with  the  camera  at  a  given 
observation  site  points  in  the  direction  of  the  optic  axis  of  the  camera  lenses. 
The  xy  plane  is  perpendicular  to  this  axis  and  passes  through  the  image  nodal 
point  of  the  optical  system.  The  selection  of  x  and  y  directions  in  this  plane  is 
arbitrary.  It  is  convenient  to  take  the  x  axis  parallel  to  one  of  the  edges  of  the 
picture  frame,  or  along  the  direction  determined  by  appropriate  fiducial  marks 
that  appear  on  each  exposure.  The  y  axis  is  chosen  to  complete  a  righthanded 
orthogonal  system. 
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As  indicated  in  die  footnote  in  the  introduction,  it  will  b-3  assumed  that  an 
ideal  camera  system  if  used  at  each  observation  site.  This  means,  for  instance, 
that  the  film  plane  is  strictly  parallel  to  the  xy  coordinate  plane  and  is  normal  to 
the  optic  axis.  It  also  presupposes  that  the  optical  system  is  totally  free  of  dis¬ 
tortion,  field  curvature.,  and  chromatic  aberration.  When  these  conditions  are 
satisfied,  distance  measurements  on  the  image  plane  can  be  readily  converted 
into  angular  measurements  relative  to  the  topocentric  coordinate  system  defined 
above.  From  Figure  3  it  follows  that : 
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Figuie  3.  Camera  Coordinate  System 
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The  unit  vector  c,  along  the  line  of  sight  to  a  point  whose  image  is  found  at  P’,  is 
therefore 


A  A  .  .  A  .  A 

c  =  x  (sin a  cosn)  +  y  (sin  a  sinn)  +  z  cos  a  . 


(8) 


The  coordinates  a  and  rj  specify  a  direction  in  space  with  respect  to  a  camera 
system  defined  by  a  camera  whose  effective  focal  length  is  f.  A  space  point  r 
units  distant  from  the  origin,  in  the  direction  of  the  vector  c,  is  represented  by 
the  position  vector 


-*•  A 

r  *  rc 


(9) 


This  completes  oar  discussion  of  the  basic  coordinate  systems  needed  to 
describe  the  positions  of  points  on  the  clouds  and  trails  used  for  the  study  of  upper 
atmosphere  phenomena.  The  auxiliary  systems  mentioned  earlier  will  be  intro¬ 
duced  in  Sec.  7.  We  will  now  relate  these  systems  in  order  to  refer  all  observa¬ 
tions  to  a  single  system  before  we  perform  the  data  manipulation  that  will  enable 
us  to  extract  physically  significant  information  from  the  raw  photographic  data. 


3.  COORDINATE  TRANSFORMATIONS 

The  position  of  a  point  in  space  is  specified  by  coordinates  referred  to  a 
system  that  is  to  a  great  extent  selected  arbitrarily  but  is  often  governed  by 
tradition  or  preferred  for  reasons  of  convenience  or  simplicity.  Since  photo¬ 
graphic  data  to  evaluate  the  motion  and  growth  of  chemical  releases  in  the  upper 
atmosphere  are  taken  from  stations  that  move  with  the  earth,  it  has  been  found 
simpler  and  more  convenient  to  refer  the  coordinates  of  points  on  clouds  and 
trails  to  the  geocentric  system  defined  in  Sec.  2.  To  accomplish  this  it  is 
necessary  to  relate  the  geocentric  and  the  camera  systems,  that  is,  we  must  know 
the  geocentric  orientation  of  the  cameras.  An  accurate  and  reliable  way  of  deter¬ 
mining  this  orientation  consists  in  photographing  the  stellar  background  shortly 
before  a  chemical  release  takes  place.  In  view  of  the  fact  that  the  star  coordinates 
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generally  listed  in  star  catalogs  ire  right  ascension  and  declination,  it  also 
becomes  necessary  to  find  the  transformation  from  equatorial  to  geocentric 
coordinates.  This  transformation  is  a  simple  two-dimensional  rotation  because 
the  two  systems  have  a  common  origin  and  coincident  z  axes.  The  transformation 
from  camera  to  equatorial  or  geocentric  coordinates  is  more  involved  and  in  its 
general  form  can  be  decomposed  into  a  three-dimensional  rotation  followed  by  a 
translation. 

A  simplification  occurs  in  photography  of  the  stellar  background.  The  very 
great  distances  involved  validate  the  assumption  that  stars  are  essentially  at 
infinity,  and  the  error  introduced  by  a  parallel  translation  of  the  camera  system 
from  the  earth's  surface  until  its  origin  coincides  with  the  center  of  the  earth  is 
therefore  negligible.  The  transformation  between  the  camera  and  equatorial  or 
geocentric  coordinates  consequently  becomes  a  general  rotation  of  axes. 

It  is  well  known  that  rotations  can  be  regarded  as  transformations  effected  by 
operators  that  can  be  represented  by  3x3  orthogonal  matrices.  The  general 
matrix  can  be  generated  by  multiplication,  in  an  appropriate  order,  of  matrices 
representing  simple  rotations  about  a  coordinate  axis.  These  rotation  matrices, 
denoted  by  R.(0)— the  subscript  refers  to  the  axis  of  rotation,  with  i  =  1,  2,  3 
denoting  the  x,  y,  z  axes,  respectively— have  elements  a„  that  satisfy  the  condi¬ 
tions  (Kaula,  1966): 


j  -  1  =  i  (mod  3)  ; 
k-  2  =  i  (mod 3)  ; 


a..  =  1,  a..  =  a. .  =  cos  0  ; 

n  *  )]  kk 


a..  =  a,.  =  a..  =  a.  .  *  0,  a.,  =  -a,  .  =  sin 0 

ij  ji  ik  ki  jk  kj 


(10) 


The  signs  apply  to  righthanded  systems,  for  counterclockwise  rotations  as 
viewed  from  a  point  on  the  positive  half  of  the  axis  looking  toward  the  origin.  For 
example,  for  rotations  about  the  z  axis,  we  have 


/  cos  0 

sin© 

0  \ 

Rz<0>  =  | 

-sin  0 

\  o 

cos  0 

0 

0 

1  ) 

\ 

(ID 

In  the  particular  case  of  conversion  from  equatorial  to  geocentric  coordinates 
we  have  angle  0  =  GHA  T.  The  connections  between  the  various  angles  are  shown 
in  Figure  4,  which  depicts  the  apparent  motion  of  the  equatorial  frame  as  viewed 


H 


[  \ 

V  w 


(GREENWICH) 


Figure  4.  Relation  Between  Geocentric  and  Equatorial  Systems 


by  an  observer  in  the  geocentric  system.  It  is  conventional  to  measure  hour  angles 
positive  westward,  from  the  local  meridian  to  the  hour  circle  of  the  star.  The 
hour  angle  of  the  equinox  is,  by  definition,  the  sidereal  time.  Let  h  and  hG  be  the 
local  and  Greenwich  hour  angles  of  a  star  whose  right  ascension  is  a.  Similarly, 
let  t  and  be  the  corresponding  sidereal  times.  Then,  from  Figure  4,  we  have 
that  at  the  instant  a  star  is  observed, 

r  =  LST  =  h  +  a  ,  (12) 

=  GST  =  h„  +  a  =  r  +  X.  (13) 
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and 

hG  =  t  +  X  -  a  .  (14) 

It  is  clear  from  the  diagram  that  we  must  set  6  =  Tq. 

The  transformation  connecting  the  camera  and  the  equatorial  systems  will  now 
be  discussed.  As  indicated  earlier,  for  the  stellar  background  we  can  use  a  pure 
rotation.  Figure  5  shows  (a)  the  basic  geometry  in  these  circumstances,  and 
(b)  how  the  general  rotation  can  be  decomposed  into  three  simple  rotations: 

(1)  about  the  z  axis  (earth's  polar  axis),  (2)  about  the  y'  axis,  and  (3)  about  the 
z"  axis.  The  rotation  angles  are  therefore  the  well-known  Euler  angles. 


Figure  5.  Relation  Between  Camera  and  Equatorial  Systems 
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Let  or^  and  6c  be  the  right  ascension  and  declination,  respectively,  of  the  point 
on  the  celestial  sphere  that  marks  the  direction  in  space  along  the  camera  optic 
axis.  Let  a,  rj'»  and  A  be  one  side  and  two  angles  of  the  spherical  triangle  formed 
by  the  camera  axis  intercept  C,  the  star  S,  and  the  pole  P.  On  the  photographic 
plate  the  star  image  appears  as  in  Figure  5(c).  If  the  star  coordinates  (right 
ascension  and  declination)  are  a  and  6,  respectively,  we  have: 


cost)  sin  a 


sinrj  sino  J  =  "  5  c)  ^z(ac)  I  s*norcos6  j» 


cos  a  cos  6 


where  Rz(6e)  and  Rz(ac)  are  given  by  Eq.  (11),  with  0  =  fie  and  0  =  c*c,  respectively, 
and 


(16) 


The  column  vectors  will  be  recognized  as  alternative  representations  of  the  unit 
vector  along  the  direction  of  the  star  in  the  two  coordinate  systems.  The  compo¬ 
nents  were  displayed  in  Eqs.  (8)  and  (1).  Equation  (15)  can  be  simplified  for 
fie  =  0  or  fie  =  ir/ 2.  In  the  former  case  it  can  be  written  in  expanded  form  as 


cos  6  cos  a  cos  or  sinfi  +  cos  5  sin  a  sin  a  sin  6  -sin  6  cos  6 
c  c  c  c  c 

=  cos  77  sin  a  , 


(17a) 


-cos  fi  cos  a  sin  a  +  cos  6  sin  a  cos  a  =  sin  n  sin  a, 
c  c 


(17b) 


cos  6  cos  a  cos  «c  cos  6c  +  cos  6  sin  a  sin  a c  cos  6C  +  sin  fi  sin  6c  -  cos  a  ,  (17c) 

which  can  be  solved  explicitly  for  fi  and  a  in  terms  of  r/  and  a.  We  can  easily 
verify  the  results 

.  ,  ,  _ sin  r]  sin  o _ 

n  to  ac>  cos cr  cos 6c  +  sino  +  sin 6c  cost) 


(18a) 
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and 


sinfi 


sin  6  cos  a  -  cos  6  sin  a  cos  q 
c  c 


(18b) 


These  expressions,  or  the  more  general  ones  that  can  be  derived  from  Eq.  (15) 
with  6e  /  0,  allow  us  to  compute  the  right  ascension  and  declination  of  an  object 
whose  image  coordinates  have  been  measured  and  used  for  evaluating  the  angles 
oand  q  [Eqs.  (5),  (6),  and  (7)].  This  is  seldom  the  way  Eq.  (15)  is  used,  how¬ 
ever.  The  primary  reason  for  establishing  Eq.  (15)  is  to  have  the  means  of 
determining  the  orientation  of  the  camera  system  from  photographs  of  the  stellar 
background  at  the  observation  site.  Details  of  the  procedure  are  given  in  Sec.  4. 
Let  us,  for  the  moment,  assume  that  a  and  &c  are  known.  It  follows,  then,  that 
in  the  equatorial  system  a  unit  vector  n  in  the  direction  of  the  camera  axis  has 
the  form 

A  A  A  A 

(n)E  =  x£  (cos  ^coso,)  +  yE(cos  sino^)  +  zEsin6c  .  (19) 

The  same  vector  can  be  referred  to  the  geocentric  system  by  direct  application  of 
Eqs.  (12)  and  (14).  We  then  have 

(n)G  =•  Xq  |cos  6c  cos  (tg  -  ac)J  +  yG  |cos  6c  sin  (tg  -  ac)J  +  zG  sin  6c  .  (20) 


The  z  axis  of  the  camera  system  coincides  with  the  optic  axis.  The  unit  vector 
zc  is  therefore  identical  with  the  unit  vector  n,  and  its  geocentric  representation 
is  also  given  by  Eq.  (20). 

When  points  at  finite  distances  from  the  camera  are  observed,  the  general 
transformations  must  be  used  but  knowledge  of  the  orientation  of  the  z  axis  of  the 
camera  system  with  respect  to  the  geocentric  system  is  prerequisite  to  the  trans¬ 
formation.  Consider  the  situation  depicted  in  Figure  6.  The  geocentric  coordinates 
are  represented  by  nonsubscripted  variables.  Unit  vectors  associated  with  the 
camera  system  are  identified  by  the  subscript  c.  Their  components  are  to  be 
referred  to  the  geocentric  system.  Let  c  be  a  unit  vector  along  the  line  of  sight 
to  the  point  P.  If  the  distance  SP  =  a,  then  the  vector  SP  =  ac.  The  vector  c, 
referred  to  the  camera  system,  can  be  written  in  the  form 

c=pzc  +  qOc,  (21) 

where 


(22) 


L4 


-y 

Figure  6.  Relation  Between  Camera  and  Geocentric  Systems 


Clearly,  the  vector  vc  is  contained  in  the  plane  determined  by  the  line  of  sight  and 

the  camera  optic  axis  and,  because  of  Eq.  (22),  is  in  the  camera  system  xy  plane. 

Figure  7  represents  a  photograph  of  P.  The  fiducial  marks  determine  plate 

axes  from  which  the  distances  x'  and  y'  are  measured.  For  a  rotation  angle 

6e  =  0,  the  plate  axes  coincide  with  the  x  y  axes,  which  define  the  E-W  and  N-S 

directions  on  the  plate.  Clearly,  the  unit  vector  x  ,  which  points  toward  the  east, 

a  a  7. 

is  proportional  to  the  vector  product  of  z  and  zc,  that  is. 


If  for  simplicity  we  write  Eq.  (20)  in  the  form 

-A  AAA  A 

(n)G  =  zc  =  x^  +  yn2  +  zn3  , 


(24) 


/ 


/FIDUCIAL  MARKS 


Figure  7.  Coordinate  Systems  on  the  Film  Plane 


where 


n,  =  cos  6  cos  (t„  -  a  ) , 

1  C  Lr  C 


n„  =  cos6„  sin(T_,  -  a  )  , 
£,  C  la  C 


n3  =  sin6c  ,  i 


Eq.  (23)  becomes 


/  2  ,  2 
ni  +  n2 


(~™2  +  H) 


The  unit  vector  along  the  third  axis  of  the  cartesian  camera  system,  y  ,  is  simply 

i  c 

the  vector  product 
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or,  in  component  form. 


(27) 


Since  det  R  =  1  and  the  matrix  is  real,  the  inverse  transformation  is  effected  by 
the  transposed  matrix 
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/A 

R*1  =  RT  = 

-nin3/A 

Vnl 

n2 

n3  / 

(31) 


The  matrix  (30)  can  be  obtained  directly  from  Eq.  (15)  with  fie  =  0,  and  from 
Eq.  (11)  with  6  -  r Evidently, 


EsVtg,r2I(Vr;1<|-  *c» 


R2(TG  "  "c1  Rj'  'f  '  V  • 


(32) 


The  vector  c,  as  given  by  Eq.  (29),  is  obtained  by  performing,  in  addition  to  the 
rotation  R,  a  rotation  about  the  zc  axis  through  the  angle  -Se.  In  other  words. 


/s\ 

* 

y 

VJ 


<'c  -  «C> 


!  cos  n  sinu^ 
sin  t)  sin  a 

coso  J 


(33) 


The  vectorial  derivation  previously  given  serves  to  identify  x£  and  yc  as 
vectors  that  point  east  and  north,  respectively,  but  for  routine  computations  with 
high-speed  digital  equipment,  Eq.  (33)  is  to  be  preferred.  This  is  because  it  can 
be  implemented  by  means  of  subroutines  of  general  applicability  in  the  triangulation 
procedure. 

Having  referred  the  unit  vector  along  the  line  of  sight  to  the  geocentric  system 
we  are  now  able  to  write  down  the  expression  for  the  position  vector  OP  =  p  . 
Clearly, 

p  =  r  +  ac  ,  (34) 

where,  we  recall,  r  is  tne  vector  that  fixes  the  position  of  the  origin  of  the  camera 
system  at  site  S.  It  can  be  recognized  that  the  addition  of  r  to  the  vector  ac  cor¬ 
responds  to  the  translational  part  of  the  transformation  that  relates  position 
measurements  in  the  camera  system  to  position  measurements  in  the  geocentric 
system.  Later  on  we  shall  see  how  r  is  determined  from  geodetic  data  appropriate 
to  the  observation  site.  For  the  moment  it  suffices  to  assume  that  the  vector  r  is 


specified  by  its  cartesian  components 


-♦  A  A  A  A 

r  =  rr  =  xx  +  yy  +  zz 

=  r  (x  cos  +  y  cos  y.^  +  z  cos  y g)  .  (35) 


The  direction  cosines  cosy,  can  be  expressed  as  functions  of  the  geocentric 
latitude  9'  and  the  longitude  X .  The  following  can  easily  be  verified: 


/  cos  yj  \ 


cos  y„ 


\cosy3  ) 


f'\ 


■B,1  (-»)»,'  <2  - 


W 


/cos  X  COS  <t>'\ 


sin  X  cos 


V 


sin< 


(36) 


Equation  (36)  is  consistent  with  the  convention  that  western  longitudes  and  southern 
latitudes  are  negative. 

Equations  (20),  (29),  and  (34)  solve  the  problem  of  expressing  in  terms  of 
geocentric  coordinates  the  position  of  a  point  at  a  finite  distance  from  a  triangula¬ 
tion  site  whose  position  and  orientation  are  known.  In  Sec.  4  we  discuss  the  means 
to  determine  that  orientation. 


4.  ORIENTATION  OK  THE  CAMERA  SYSTEM 

Equation  (15)  relates  the  right  ascension  and  declination  of  a  star  to  the 
angular  coordinates  a  and  r,  that  correspond  to  the  orientation,  in  the  camera 
system,  of  the  line  of  sight  to  the  star.  The  transformation  (15)  depends  on 
three  parameters,  of  which  «c  and  6c  are  identified  as  the  right  ascension  and 
declination  of  the  point  on  the  celestial  sphere  toward  which  the  camera  optic 
axis  is  directed,  and  6e  corresponds  to  a  rotation  of  the  film  plane  about  the  optic 
axis. 

We  will  now  set  up  a  procedure  to  determine  a c,  6  ,  and  6e  from  photographic 
records  of  stars.  Properly  identified  by  comparison  with  charts  from  a  star 
atlas,  they  will  supply  the  data  from  which  the  column  vectors  on  the  righthand 
side  of  Eq.  (15)  can  be  constructed.  Measurements  of  the  corresponding  star 
images  on  the  film  will  provide  the  data  that,  coupled  with  Eqs.  (5),  (6),  and  (7), 
will  allow  evaluation  of  the  components  of  the  column  vectors  on  the  lefthand  side 
of  Eq.  (6).  Since  the  position  of  each  star  can  be  described  by  only  two  independent 
equations  involving  its  right  ascension  and  declination,  identification  of  a  single 
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cos  6  cos  6  cos  (a  -a)  +  sin  6  sin6=cosn.  (39c) 

c  c  c  — 

Equations  (39)  can  be  written  in  the  equivalent  form 

cos  6  sin  ( «c~a )  +  sin  a  sin  (q  +  fie)  =  0  ,  (40a) 

sin6c  cos  6  cos  (<*c~a)  ~  cosfic  sin  6  =  sina  cos  (q  +  6e)  ,  (40b) 

cos6c  cos  6  cos  {ac~a)  +  sinfic  sin  6  =  cos  a  .  (40£) 

The  added  rotation  of  the  film  plane  about  the  optic  axis  appears  as  a  correction  in 
the  angle  q,  which  measures  the  azimuth  of  the  star  image  from  the  reference 
direction  determined  by  the  fiducial  marks. 

Equations  (40b)  and  (40c)  can  be  used  to  express  sinfic  and  cos  6c  as  functions 
of  «c  and  6e.  It  can  readily  be  seen  that 
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8in6c  =  2^  [ 

sin 2 6  cos  fac~o )  +  sin2o  cos  (q  +  6e)  J 

*tl 

sin 6  cos  a  +  cos  6  sina  cos  (o+  6e)  cos  fa  -a)] . 

(4  la) 

cos6c=2i;  | 

cos  26  +  cos  2a  J  =  j^cos  6  cos  a  cos  faQ  -  a)  - 

-  sin 6  sina  cos  (rj  +  6e)  J  , 

(41b) 

where 

Aj  =  coso  cos 6  cos  fac~a)  +  sino  sin 6  cos  (n+  6e)  , 

2 

c  l-cos  6 

sin2  fac  -  a)  . 

(41c) 

Consequently, 


sin  26  cos  fa  ~a)  +  sin2o  cos  (rj  +  6e) 
ten6c  cos  26  + cos  2o 

sin  6  coso  +  cos  6  sin  a  cos  (rj  +  6e)  cos  fa  -a) 

—  —  ■■  ...  .  .  - - - -  JS —  -  -  0 

cos  6  cos  a  cos  fac  -a)  -  sin  6  sin  a  cos  (n  +  6e) 
We  now  consider  Eq.  (40a),  rewritten  in  the  form 


sin<?c  cos  6  cos a  -  cosac  cos  6  sina  =  -sino  sin  (rj  +  6e)  .  (43) 

Then  we  evaluate  Eq.  (43)  for  each  of  two  stars  and  solve  for  sino,  and  cos«c> 
With  the  subscripts  1  and  2  corresponding  to  the  first  and  second  stars,  respect¬ 
ively,  we  have 

(44a) 

(44b) 

where 


sin  a 


sinu2  sinafj  -  sinUj  sinOg 

Sin(ar1  -  ar2) 


smu  cos  a,  -  sinu,  cosa, 

cose  = - 2  ,  1 - - - 

c  sin  (a1  -  a2) 


aino.  \ 

-^T.  j  Sinope). 

i  ' 


sin  u^  = 


(44c) 


2  2 

By  virtue  of  the  identity  sin  crc  +  cos  <*c  =  1,  we  can  write 
2  2  2 

sin  ((*1  -  a 2)  *  sin  +  sin  -  2  sinu^  sinu^  cos  (ar^-crg) 


2  1  /sinai  \2  /  \ 

3in  u|  =  2~  \cos6.  /  \  1  "  cos2,Ji  cos26e  +  s'n2rji  sin26c  j  , 

1  /8inai  \/sina2\r 

sinuj  sin*,  =  j  -  cos  26c  cos  (rjj  +  n2>  H 

+  stn26c  sinfrjj  +  rj2> ] 


Consequently,  we  can  write 


Sj  cos  26c  -  S2  sin  26c  =  Sq  , 


where 


2 

/sina,  \  /  sin  a,  \  /  sina  \ 

si =  cos2TJi  v^3i Tj )  •  c°8(arra2 


/sin a„  \  ^ 

+  ( - j-  )  cos  2q„  , 

\cos62  /  '2 


2 

/sina,  \  /sinajWsinag  \ 

S  =  sin2rj  ( - r~  )  -  2  (  - )( — — t  )  sin  (rj, +rj„)  cos  (ar, -ar„)  + 

2  lycose,/  \cob  b  ^)\cob  b ^  )  1  2  1  2 

/sina  \2 

c°s2i72  , 

/sina1\2  /sina,  \  /  sin a_  \  /sina  \2 

So  \<ms6^)  "  2(cos 6,  ) \cosl^  )  COS  (®1  “V  +  \cos“^  ) 


o 

-  2  sin  (a-j  -  a2)  . 


From  Eq.  (45)  it  follows  that 


cos  (26c  +  tp)  =' 


Sl  +  S2  ’ 
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with 


ip  s  tan 


From  this  equation  we  find  the  following  expressions  for  6c: 


From  Eqs.  (46)  it  is  apparent  that  the  rotation  angle  6c  vanishes  for  =  Sq. 
This  leads  to  the  condition 


( 


sinOj  sinrjj 
cos  6^ 


2  ni+ilo  "l  2 

X  sin  — ^ - sinrjj  sinrjg  +  sin  (afj  - a 2)  , 

which  can  be  used  in  conjunction  with  Eqs.  (37)  and  (38)  to  check  the  consistency 

of  the  data  from  pairs  of  stars  for  known  6c  =  0. 

Having  determined  the  angle  6c,  which  measures  the  rotation  of  the  xy  plane 

of  the  camera  system  from  the  N-S  direction,  we  can  use  Eqs.  (44)  and  (42)  to 

find  the  angles  a  and  6  that  identify  the  orientation  of  the  optic  axis  (or,  equi- 
c  c 

valently,  the  orientation  of  the  camera  system  z  axis).  We  have 


(47) 


It  should  be  remarked  that  the  statements  made  at  the  outset  concerning  the  need 
to  use  data  from  at  least  a  pair  of  stars  in  order  to  evaluate  the  angles  <*c,  6C,  and 
fie  is  fully  consistent  with  the  system  of  Eqs.  (40).  We  have  here  three  equations, 
which  in  ordinary  circumstances  should  be  sufficient  to  solve  for  the  three  un¬ 
knowns.  Closer  examination,  however,  discloses  that  they  are  not  independent. 
Unless  previous  information  is  available  about  the  value  of  any  one  of  the  unknowns, 
a  sufficient  number  of  independent  equations  can  only  be  obtained  by  adjoining  to  the 
set  of  equations  from  one  star  the  corresponding  set  for  a  second  star. 


5.  TRIANGULATION  PROCEDURE  FOR  POINT  RELEASES 

Consider  now  the  problem  of  determining  the  geocentric  coordinates  of  the 
point  in  space  that  marks  the  position  of  a  chemical  puff.  To  locate  the  puff, 
photographs  are  taken  from  two  or  more  stations  whose  geodetic  coordinates  are 
known.  Let  us  assume  that  the  geocentric  coordinates  of  the  stations  have  been 
found  by  an  appropriate  transformation,  *  Construct  the  position  vectors  r.  that 
fix  the  location  of  the  stations  in  the  geocentric  system.  They  have  the  general 
form  specified  by  Eqs.  (3)  and  (4),  that  is, 

*  This  transformation  is  discussed  in  Sec.  7. 
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r.  = 
1 


A 

r.g. 


where 


g^  =  x  (cos  X  cos  tfj)  +  y  (sin  V  cos  $j)  +  z  sin  , 

and  $!,  X  denote  the  geocentric  latitude  and  longitude  of  the  i_th  station.  At  each 
observation  site  a  camera  is  mounted  with  its  axis  oriented  in  a  direction  specified 
by  the  unit  vector  n,  whose  geocentric  components  were  determined  in  Sec.  3, 

Eq.  (20),  such  that 


n  *  x  cos  fi  cos  (t„  -  a  )  +  y  cos  6  sin  (t „-a  )  +  z  sin  6  . 

c  G  c  c  G  c  c 

There  will  be  one  such  n  vector  for  each  coordinate  system  associated  with  a 
camera.  To  avoid  an  overabundance  of  subscripts  this  equation  shows  the  general 
form  of  n.  When  it  becomes  necessary  to  identify  a  particular  n,  we  will  attach 
the  subscript  and  write  n..  The  expanded  form  for  n.  is  computed  from  the  or  , 

6c,  and  rG  appropriate  to  the  camera  and  observation  site  under  discussion. 

Let  us, now  consider  the  situation  depicted  in  Figure  8.  A  point  release  at  P  is 
photographed  from  sites  Sj  and  S2-  It  was  shown  in  Sec.  3  that  along  the  line  of 
sight,  a  unit  vector  c  is  expressed  in  geocentric  coordinates  by  Eq.  (29),  repeated 
here  for  convenience: 


c  =  x  n^coso  -  °  |\i2  cos  (b  +  6e)  +  njn3  sin  (rj+  6e)j  j  + 


.  a  )  ,  sino 

+  y  (  n„  cos  a  + 


^nj  cos  (rj  +  fie)  -  n^g  sin  (q+  6e)J  j  + 


+  z  £  ng  cos  o+6  sin  or  sin  (n  +  fie)  J 


[(29)] 


where 

.  /  2.2 
A  =Vnl+n2  • 

Let  us  introduce  the  notation: 


OP  =  p 
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Figure  8.  Triangulation  on  Point  Releases 


The  triangulation  problem  consists  of  finding  an  expression  for  the  vector  p.  From 
Figure  8  it  is  evident  that  we  can  write 


-V  A 

r1  +acx  , 

(50a) 

■+  .  A 

r2  bc2  * 

(50b) 

26 


By  a  scalar  multiplication  of  Eqs.  (50a)  and  (50b)  we  obtain  for  the  magnitude  of 
p  the  expression 


/-*■  A  A  1  a  A  \  l/2  .... 

p  =  (  ri»r2  +  ar2*ci  +  bri*c2  +  ab  C1*C2  )  '  *  (51> 

Vector  multiplication  leads  to  the  relation: 

0  =  pxp  =  +  a(CjXr2)  +  b^xc^  +  abfCjXCg)  .  (52) 

Take  now  the  scalar  product  of  Eq.  (52)  and  or  c^.  This  procedure  leads  to  the 
following  explicit  expressions  for  the  scale  factors  a  and  b: 


a  -*  -* 

Vrl*r2 

A  A 

r2'clXc2 


-*■  A 

IVr2XC2 

A  *>  A 

cl*r2XC2 


(5?a) 


A  — , 

-  VriXr 


A  A 

rl*ClXc2 


2  _ 


-*■  w  A 

VriXci 

A  -»  -,A 

C2*rlXC! 


(53b) 


The  consistency  of  the  scheme  can  be  checked  by  virtue  of  the  relation 

,2  _  ,.a  a  .2 _ 2  ,  .2 


I*d  =  d/  =  (bc2  -  ac =  3l  +  D  -  2ab  (Cj'Cg) » 


(54) 


2 

in  which  d  represents  the  square  of  the  distance  separating  the  two  observation 
sites.  This  distance  is  often  directly  available,  having  been  ascertained  when  the 
region  where  the  observation  sites  are  located  was  surveyed.  Alternative  expres¬ 
sions  for  the  position  vector  p  are  obtained  when  values  of  a  and  b  are  substituted 
in  Eqs.  (50a)  and  (50b).  Since  there  is  no  a  priori  reason  to  favor  one  observation 

-V  I 

site  over  the  other,  it  is  advisable  to  define  p  as  a  symmetric  combination  of 
Eqs.  (50a)  and  (50b),  that  is, 

P  2  (rl+r2+aCl  +  bC2)  '  <5S> 

In  combination  with  Eq.  (51)  this  expression  can  be  used  to  define  a  unit  vector  p 

i 

in  the  direction  of  OP,  the  position  vector  of  the  release.  Knowledge  of  this  unit 
vector  allows  immediate  calculation  of  the  geocentric  latitude  of  the  release. 
Clearly, 


p»z  =  cos(t:  ~  4>'  )  sin  <j>} 

'2  p  /  p 


(56) 
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where  4^  denotes  the  geocentric  latitude  of  P.  Similarly,  the  longitude  of  the 
release  point  can  be  computed  from  the  expression 


tan  A  *  ~~X-  • 

PA  A 
Cx  a  Y 


If  only  the  cartesian  geocentric  coordinates  of  the  point  P  are  required, 

Eq.  (55)  is  all  that  is  needed.  Equations  (51),  (56),  and  (57)  allow  parametrisation 
of  these  components,  thus  serving  to  locate  P  in  the  geocentric  system  by  means 
of  spherical  coordinates.  This  is  particularly  useful  because  geodetic  coordinates 
can  be  obtained  from  thic  form  very  readily.  (  The  transformation  formulas  will 
be  discussed  in  Sec.  7.  ] 

The  position  vector  p,  which  locates  the  point  P  in  the  geocentric  frame  of 
reference,  is  uniquely  determined  when  independent  observations  from  two  or 
more  stations  are  available.  Equations  (55)  and  (53a,  b)  solve  the  two-station 
problem.  It  can  easily  be  shown  that  for  N  stations  the  position  vector  p  is  given 
by  the  expression: 


Z(V*A)- 


where 


r j  =  position  vector  of  jth  camera  site 

Cj  =  unit  vector  along  the  direction  of  the  jth  camera  optic  axis 

=  scale  factors;  for  the  jth  line  of  sight,  given  by  the  vector  relations 


V 


A 

ri*ci 

i=l _ 

N 

£  ?ixSi 


There  is  some  arbitrariness  in  the  construction  of  a  check  equation  similar 
to  Equation  (54)  for  the  two-station  case.  A  symmetric  form  can  be  built  by  taking 
the  stations  in  a  close  cycle  and  computing 
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(60) 


where  the  index  N+l  is  to  be  identified  with  the  index  1,  a  condition  necessary  to 
complete  an  arbitrary  closed  cycle  of  stations  with  no  cross  links. 

Before  concluding  this  section  a  reminder  is  in  order  about  the  significance  of 
Eqs.  (58)  and  (59).  It  has  been  implicit  that  the  analysis  that  led  to  Eq.  (59) 
assumes  that  ail  lines  of  sight  intersect  at  the  point  P,  which  coincides  with  some 
structural  feature  of  the  release  uniquely  identifiable  from  all  observation  sites. 

Then,  Eq.  (58)  represents  the  position  vector  obtained  by  making  use  of  all 
observations  symmetrically,  with  equal  weights.  In  a  sense,  it  is  a  mean  value 
of  equally  weighted  observations.  The  more  interesting,  but  more  involved 
problem  of  deriving  expressions  for  the  probable  error  and  standard  deviation 
associated  with  the  mean  value  given  by  Eq.  (58)  will  be  the  subject  of  a  separate 
report.  Thus,  the  simplicity  of  our  approach  is  revealed  at  all  stages,  unencumbered 
by  the  complexities  of  a  full  statistical  analysis. 


6.  TR1ANGLLAT10N  ON  CONTINUOUS  TRAILS 

For  the  triangulation  scheme  that  was  discussed  in  Sec.  5  it  is  assumed  that 
the  images  of  a  distinct  cloud  feature  can  be  individually  identified  on  photographs 
from  different  sites.  This  is  not  possible  when  the  chemical  release  takes  the 
form  of  a  continuous  trail.  For  example,  from  a  pair  of  observation  sites  the  two 
trail  images  have  the  general  appearance  shown  in  Figure  9.  This  figure  clearly 
indicates  that,  except  possibly  for  the  trail  end,  it  is  not  usually  possible  to 
establish  a  unique  correspondence  between  a  given  point  Pj  on  trail  image  1  and 
some  point  P2  on  trail  image  2. 

To  match  point  pairs  on  the  two  trails  it  is  necessary  to  use  a  more  general 
computational  scheme  for  the  triangulation,  one  that  allows  the  position  of  point  Pg 
along  the  second  trail  image  to  be  varied  until  some  criterion  for  optimum  match 
is  satisfied.  This  requires  a  procedure  for  moving  along  the  trail  image  and 
scanning  over  short  sections  of  it  to  locate  the  optimum  pairing.  A  very  convenient 
way  to  accomplish  this  is  to  replace  the  actual  trail  image  by  a  two-dimensional 
curve,  which  is  represented  analytically  by  a  parametric  cubic  spline  (Ahlberg, 
Nilson,  and  Walsh,  1967).  For  computer  implementation,  the  scan  can  be  effected 
by  systematically  incrementing  the  spline  parameter.  The  merit  function  that 
measures  the  quality  of  the  match  is  computed  after  each  increment.  In  general, 
it  will  decrease  as  the  scan  proceeds,  reach  a  minimum,  and  then  increase 
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Figure  9.  Luminescent  Trail  as  Observed  From  Two  Different 
Triangulation  Sites 


(Bekey  and  McGhee,  1964).  If  the  increments  are  small  enough,  or  if  suitable 
interpolation  procedures  are  used,  the  minimum  value  can  be  determined  very 
precisely. 

It  is  not  the  purpose  of  this  report  to  dwell  on  this  aspect  of  the  triangulation 
problem.  Rather,  let  us  try  to  establish  the  analytic  form  of  the  merit  function 
and  assume  that  adequate  procedures  are  available  for  an  effective  determination 
of  the  optimum  match. 

Consider,  then,  the  situation  depicted  in  Figure  10.  The  vector  diagram 
shows  the  essential  geometry  involved.  Clearly, 


t  = 


A 

=  dd 


and 


(61) 


A 

ac 


1 


-  bc£  +  del  =  0  . 


s 


(62) 
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To  determine  the  scale  factors  a  and  b  we  define  s  *  |  s  |  as  the  minimum  dis¬ 
tance  between  the  lines  of  sight.  We  can  thee  write 


-*■  ,  .  A  A 

s  *  dd  +  aCj  -  bc9  , 

(63a) 

+  A  A 

6«Cj  *  S.C2  *  0  , 

(63b) 

whence  we  obtain  the  following  expressions  for  a  and  b 


[  d*c  -  (3.C-)  (C-C  )  1 

'H'  i  +  (vV2  -I’ 


(64) 


=  d 


rA  A  A  A  A  A  , 

d  *  c2  -  ^’Ci>  <c1‘c2)l 

1-  1  +  (c^Cg)2  J 


(65) 


31 


A  vector  p  can  now  be  constructed  from  the  origin  to  the  midpoint  of  the  line  seg¬ 
ment  PjP2*  virtue  of  conditions  (63b),  PjP2  is  shortest  distance  between 
the  lines  of  sight  and,  clearly. 


■*  ■*  .  A  1  -» 

p  =  ri  +aci  -JB  , 


4  ,  A  1-4 

r2  +bc2  +28 


Consequently, 

^*2(?l+?2+a^l+b^2)  •  (66) 

Equation  (66)  is  identical  with  Eq.  (55),  which  was  obtained  for  point  releases.  The 
disappearance  of  any  reference  to  s  in  Eq.  (66)  indicates  that  the  role  of  s  is  to 
supply  a  quantitative  measure  of  the  closure  error  incurred  when  the  lines  of  sight 
and  the  base  line  connecting  the  observation  sites  are  taken  as  the  sides  of  a 
triangle  with  one  vertex  marking  the  release  point.  For  point  releases,  any  dis¬ 
crepancy  between  the  left-  and  righthand  sides  of  Eq.  (54)  serves  to  characterize 
this  closure  error.  By  modifying  the  procedure  used  to  evaluate  a  and  b,  it  can 
easily  be  verified  that  Eqs.  (64)  and  (65)  can  be  obtained  for  point  releases.  In 
fact,  these  equations  emerge  instead  of  Eqs.  (53a)  and  (53b)  if  the  evaluation  of 
a  and  b  is  based  on  the  triangle  SjSgP. 

The  closure  error  reduces  to  zero  whenever  the  lines  of  sight  intersect;  how¬ 
ever,  the  occurrence  of  intersecting  lines  of  sight  does  not  mean  intersection  at 
the  point  where  the  release  has  taken  place.  To  estimate  how  far  away  from  the 
actual  release  point  the  intersection  occurs  it  is  necessary  to  have  data  from  more 
than  two  observation  sites.  The  statistics  of  the  p  vectors  will  be  discussed  in  a 
separate  report. 

Within  the  limitations  imposed  by  triangulation  from  two  sites,  we  now  construct 
the  function  that  allows  matching  points  on  pairs  of  trail  images.  Let  us  define  the 
unit  vectors  v^  and  v2  as  follows: 


a  Ci 

V1  *  sinJ3j  ’ 

A  A 

a  _dXC2 

v2  ”  sin^2 


(67a) 


(67b) 
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where  Vj  is  a  vector  normal  to  the  plane  determined  by  the  lines  SjS2  and  S^P^, 
v2  is  normal  to  the  plane  determined  by  lines  and  S,,P2,  and  sin^  is  the 
magnitude  of  the  vector  product  d  x  c  ..  It  is  evident  that  if  the  two  points  Pj  and 
P2  are  images  of  the  same  trail  point,  the  two  planes  coalesce  and  the  vectors 
Vj  and  v2  are  parallel.  When  this  occurs,  1  -  Vj«v2  vanishes.  We  can  therefore 
make  the  function 


A  A  A  A  \  / j  A  v 
c^c,  -  (d-Cj)  (d*c2) 

V[l  -  (3-c1)2]  [l  -  (d-c2)2] 


(68) 


the  merit  function  that  measures  how  close  to  coincidence  the  planes  S^SgP^  and 
SiS2P2  come.  Ideally,  this  function  vanishes  when  P2  becomes  the  second-site 
image  of  P^.  In  practice,  small  residuals  will  almost  always  be  present,  and  the 
criterion  for  a  second-site  match  must  be  the  occurrence  of  a  minimum  of 
f  [ci*  c2(t)]  as  the  parameter  t  is  varied.  This  parameter  determines  the  position 
of  the  point  P2  along  the  trail  image. 

It  is  clearly  advantageous  to  have  an  analytic  expression  to  describe  the  trail. 

In  a  perfectly  still  atmosphere,  the  trail  centerline  will  coincide  with  the  trajectory 
of  the  rocket  that  effects  the  release.  In  principle,  it  should  be  possible  to  write 
analytic  expressions  for  this  trajectory  and,  by  suitable  projective  transformations, 
arrive  at  an  equation  for  the  trail  centerline  on  the  film  plane.  Besides  the  algebraic 
complexity  that  can  be  expected  in  all  but  a  few  special  cases,  there  are  considera¬ 
tions  of  winds,  windshears,  and  turbulent  transport  processes  that  distort  the  trail 
and  often  create  sharp  twisting  and  self-intersections,  the  net  result  being  that 
the  analytic  approach  suggested  by  the  static  atmosphere  model  loses  all  sig¬ 
nificance. 

The  obvious  way  out  of  this  difficulty  is  to  employ  a  piecewise  analytic 
representation  of  the  trail  centerline  image.  Since  we  have  selected  a  parametric 
cubic  spline  representation,  the  ordinary  treatment,  which  requires  that  the  depend¬ 
ent  variable  be  a  single-valued  function  of  the  independent  variable,  has  to  be  modi¬ 
fied  to  allow  for  the  multivaluedness  introduced  by  severe  twisting  and  self- 
intersections.  An  account  of  the  parametric  spline  procedure  is  found  in  Appendix  A. 

Equations  (64),  (65),  and  (66)  solve  the  triangulation  problem  for  trails  when 
observations  are  made  from  two  stations.  As  was  the  case  for  triangulation  on 
point  releases,  it  is  possible  to  extend  the  analysis  and  develop  generalized 
expressions  for  the  N-station  case.  If  the  position  vector  p  is  computed  from  the 
equation 

N 

P  =  ^  (r-  +  a.  c.)  , 

N  Lu  j  J  J 
J=1 


(69a) 
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it  can  be  shown  that  the  scale  factors  a^  are  given  by  the  relations 


(69b) 


(70a) 


(70b) 


The  vectors  c-  are  the  unit  vectors  along  the  lines  of  sight  that  were  introduced  in 

t  A 

Sec.  3.  The  vectors  d^  are  unit  vectors  along  the  directed  line  from  the  jih  to 
the  kth  station.  The  distance  between  these  stations  is  d^.  Clearly, 


*1 


d.. 


but 


To  match  points  on  the  trail  images  taken  from  the  different  triangulation  sites 
it  is  necessary  to  construct  a  set  of  merit  functions  of  the  form  (68).  If  a  point  is 
selected  on  the  trail  image  corresponding  to  the  j_th  station,  the  spline  parameter 
values  an  the  other  images  are  determined  by  minimizing,  in  succession,  the 
functions 


1  - 


A  A  .A  0  wA  J  * 

c.  •  c.  -  (c.  •  d. .)  (c.  •  d. .) 
==j==  J_  ___i _.J__  ii 

vTi-(vV2][1_(yV2] 


(71) 


where  tj  denotes  the  parameter  for  the  jth  trail  image  when  the  point  to  be  matched 
lies  on  the  ith  trail  image.  The  set  of  parameter  values  that  is  thus  obtained  is 
used  to  determine  corresponding  points  on  the  N-l  trail  images  that  have  been 
compared  with  the_ith  image.  It  is  not  difficult  to  see  that  by  successively 
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applying  the  same  procedure  to  each  trail— taking  as  the  point  to  be  matched  the 
point  on  the  trail  determined  by  the  parameter  value  that  emerged  from  the  initial 
round— a  set  of  parameter  values  is  obtained  for  each  trail  image,  from  which  a 
mean  value  can  be  extracted,  namely, 

ft)  4  •  «« 

k«l 

These  averages  can  clearly  be  used  to  advantage  to  compute  the  line-of-sight 
unit  vectors  c^  that  appear  in  Eqs.  (69a)  and  (70).  Higher  moments  of  the  distri¬ 
butions  for  the  parameters  t.  will  be  given  elsewhere.  The  equations  developed  in 
ibis  section,  together  with  those  of  Sec.  5,  solve  the  triangulation  problems  en¬ 
countered  in  the  course  of  reducing  data  from  a  chemical  release.  They  lead  to 
expressions  for  the  geocentric  coordinates  of  the  release,  from  which  such  things 
as  wind  speeds  are  computed.  To  facilitate  comparisons  with  existing  data  it  is 
sometimes  convenient  to  express  positions  as  sets  of  coordinates  relative  to  other 
systems  of  reference.  Two  of  these  will  be  discussed  in  Sec.  7. 


7.  TRANSFORMATIONS  TO  OTHER  REFERENCE  SYSTEMS 

Although  the  equations  developed  in  Secs.  5  and  6  solve  the  problem  of 
establishing  the  position  of  chemical  releases,  these  equations  refer  coordinates 
to  a  geocentric  system.  Some  limitations  that  are  thus  imposed  can  only  be 
resolved  if  means  are  provided  to  convert  the  coordinates  to  other  reference 
systems.  Two  such  systems  (Kaula,  1966;  Mueller,  1969)  are  discussed  in  this 
section:  (1)  the  horizon  system,  which  facilitates  visualization  of  the  release 
location  with  respect  to  the  observer,  (2)  the  geodetic  system,  which  simplifies 
reference  to  geodetic  or  geographic  charts.  Since  the  transformations  to  and 
from  these  systems  can  be  found  without  difficulty  by  applying  the  appropriate 
definitions  or  the  rotation/translation  operators  introduced  in  Sec.  3,  the  trans¬ 
formations  will  be  discussed  with  a  minimum  of  detail  concerning  their  derivation. 

7.1  Geodetic  Coordinates 

Geodetic  coordinates  are  defined  with  respect  to  a  reference  ellipsoid,  a  two- 
parameter  surface  of  revolution  very  close  to  the  actual  surface  of  the  earth. 
Geometrically,  it  is  a  quadric  surface  on  which  all  curves  of  intersection  with 
planes  are  either  circles  or  ellipses.  The  ellipsoid  is  customarily  defined  by  the 
length  of  the  longest  and  shortest  axes  a  and  b,  or,  equivalently,  by  the  long 


“—to- 
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axis  a  and  one  of  the  following  auxiliary  quantities 
Flattening:  f  *  — ™- 


First  eccentricity:  e  *  \1  -  (b/a)^ 


Second  eccentricity:  e ' 


-  1 


The  ellipses  whose  generating  planes  contain  the  axis  of  rotation  are  the 
geodetic  meridians.  One  of  the  geodetic  meridians  is  selected  as  the  zero  meri¬ 
dian.  The  angle X  between  the  piane  of  the  zero  meridian  and  the  plane  of  the 
geodetic  meridian  of  a  point  P  is  the  geodetic  longitude  of  P.  Clearly,  if  the  zero 
meridian  is  contained  in  the  plane  of  the  geocentric  prime  meridian  (Greenwich), 
the  geodetic  and  geocentric  longitudes  are  identical.  The  line  PP'  (s'  Figure  11) 
perpendicular  to  the  tangent  plane  at  the  point  where  PP’  pierces  the  ellipsoid  is  a 
geodetic  normal.  The  angle  4  between  this  normal  and  a  plane  perpendicular  to  the 
axis  of  rotation  (for  example,  the  geodetic  equator)  is  the  geodetic  latitude  of  P. 
Figure  11  also  shows  the  corresponding  geocentric  latitude  4'- 


Figure  11.  Geodetic  System  of  Coordinates 
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7.3  Geodetic  to  Geoceotric  Coordiaales 


Transformations  involving  geodetic  or  horizon  coordinates  (azimuth  and 
elevation)  will  now  be  discussed.  From  Figure  11  we  can  determine  the  geocentric 
coordinates  of  a  point  whose  geodetic  latitude,  longitude,  and  elevation  above  the 
reference  ellipsoid  are  d,  X,  and  b.  respectively.  Let  the  geocentric  coordinates 
be  expressed  in  terms  of  the  parameters,  4’,  X',  and  p.  Then,  clearly, 

X'  «  X  .  (73) 

From  Figure  1 1  it  follows  that 

r  *  r  (x  cos  +  z  sin  4") ,  (74) 

h  =  h  (x  cos  $  +  z  sin  4)  .  (75) 

For  the  reference  ellipsoid  we  can  also  write 

r  *=  xa  cos  t  +  zb  sin  t  ,  (76) 


where  t  is  a  parameter. 

A  unit  vector  tangent  to  the  ellipsoid  at  M  is  given  by 
t  *  -(a  sint)x  +  (b  cost)z 

•  Va2  sin2  t+b2  cos2  t 


(77) 


But 


ii«t  ■  0  (condition  of  perpendicularity)  • 

Hence, 

tant  *  (b/a)  tan$  ,  (78) 

r  *  \/a28in2  t  +  b2  cos  2t  . 


(79) 
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Since 

tan  4"  *  (b/a)  tant  *  (b/a)2  tan  4  (80) 

and 

2  2 

(b/a)  *  1  -  e  (e  *  eccentricity)  , 
it  follows  that 


where  r  is  given  by  Eq.  (76).  From  Eqs.  (85),  (76),  and  (77),  and  the  condition  of 
perpendicularity  of  h  and  t  we  have 
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Unt  *  j,i,-g-frtt  +  ay8int 

a  p  cos  f 

a  2  1 
«  (b/a)  tan$'  +  —  co'8^i  sint .  (86) 

Since  e  «  1  and  sint  «  tant/(l  +  tan2t)^2,  the  Lagrange  expansion  leads  to  the 
expression 

tant  *  VT7U'  +  e2— + 

L  W  Vl-e2  sin2  4' 

+  i(jf  ^  — sl+Ofe6);  (87) 

2W  (1  -  e2sin2  *')2j 

Once  we  know  the  value  of  the  parameter  t,  we  can  immediately  write 

r  =  a  (r  cost  +  z  Vl  - "e2  sint)  (88a) 

i 

and 

h  =  r  (p  cos  4 1  -  a  cos  t)  +  z  (p  sin  $ '  -  a  V 1  -  e2 sin  t)  .  '  (88b) 

i 

Hence, 

i 

h=  [p2-2ap  cos  $ 'cos  t  +  Vl  -  e2  sin^1  sint  +  a2(l  -  e2sin2t)J1^2  ,  (89) 

(90) 

(91) 

The  general  transformation  from  geocentric  to  horizon  coordinates  involves  a 
translation  of  the  origin  to  the  observation  site,  followed  by  rotations  about  the  z 
and  y  axes,  respectively.  If  the  position  vector  in  the  horizon  systeih  is  p^,  it  is 
such  that  i 


,  P  sin  '  -  a  Vl  -  e2  sin  t 
=  - p  cos-T  - a  Soil - 


X  =  X' 


7.4  Geocentric  to  Horizon  Coordinates 


k 
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(■cos  E  cosA\ 

cosEsinA  ],  ■ 

sinE  / 

where  E  is  (he  elevation  and  A  is  the  azimuth  of  the  line  of  sight  to.  the  point  being 
observed.  Then 


/-cosE  cos  A\ 


PH(  cosE  sinA  J=  Ry^f  -  4'^*^) 


sinE 


/cos  cosX\ 


^cos  cos  X 
8  8 


pi  cos^'  sin*  j-r  I  cos-$'  sin* 

St  8  8 


sin 


Sin  <t>' 


I.  (92) 


where  4*  aijd  X^  &re  the  geocentric  latitude  and  longitude  of  the  observation  site 
and  r  is  the  distance  of  the  site  from  the  geocentric  origin.  It  can  be  verified 

,8  i 

that  Eq.  (92)  leads  to  the  relation 


=  1  +  ! 


cos^'  cos^'  cos  (X  -  X  )  +  sin  sin 
8  8  8 


in^'J  . 


The  term  in  the  brackets  can  be  identified  with  the  elevation  of  the  line  of 

■  .  i 

sight  to  a  target  at  a  very  great  distance  (a  star,  for  example)  and  it  is  therefore 

.  1  ■  1  . 

cohvement  to  set  > 


sinEw  *  coe4g  cos  4’  cos  (X  -  \g)  +  sin^  sin^' 


(93) 


Letting  m  ■  r  /  p  we  can  then  write 


Pit  \  2 

—  j  »  1  +  in  -  2m  sinE  , 
iP  00 


(94) 


-  m  +  sinE. 


•!sinE  = 


f  1  -  2m  sinE  >+  m 
00 


(95) 


i 
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tan  A  = 


cos  ♦'  sin  (X  -  X  ) 

co8^'  sin$'  -  sin^'  cos  6'  cos  <X  -  X  )  * 
8  S  8 


7.5  Horiui  to  Geoceairic  Cooitb  Bates  (special  case) 

The  transformation  inverse  to  Eq.  (92)  is  of  importance  only  for  very  distant 
objects  (stars).  The  transformation  from  horizon  to  geocentric  coordinates  in 
this  special  case  is  accomplished  by  application  of  the  rotations  R~Vy  -  6M  and 
R*(X  )  in  succession.  The  results  are  '  ' 


8 in  4  =  sin  6'  sin  E  +  cos  4'„  cos  E  cos  A  , 
8  8 


sinE  sinX _  cosf'  +  cos  E  (cos  X  sin  A  -  sin  X  cosAsiu^') 
X  —  8  8  8  S 

sinE  cos  X  cos  -  cosE  (sin  X  sin  A  +  cosX  cos  A  sin  6' ) 

8  8  8  S  S 


7.6  Caaera  to  Horicoa  Coordinates 

The  transformation  from  camera  to  horizon  coordinates  can  be  performed  by 
combining  a  transformation  from  a  camera  to  a  geocentric  system,  followed  by  a 
transformation  from  a  geocentric  to  a  1  zon  system.  Symbolically  we  can  write 


-cosE  cos  A' 


cosE  sinA  J  =  Ry(|-  p;)Rz^b  +  rQ  ■  c)R~j(i -  *c)  Rz('6e)l  siDrt  Bina  1*  (99*> 


'cosij  sino1 


P-  =  a. 


where  ap  is  the  distance  along  the  line  of  sight  to  the  target  being  observed.  It  is 
determined  from  an  expression  of  the  form  (59)  or  (69b). 

If  the  orientation  of  the  camera  axis  is  referred  to  the  geocentric  system,  the 
declination  6^  and  right  ascension  arn  can  be  replaced  by  an  effective  latitude 

=  fic  and  an  effective  longitude  \  s  ac  ~  tq  •  Equation  (99a)  then  take  the  form 


■cosE  cos  A' 


E  sinA  = 


cosn  sm< 


■a wppsv-r 


'-Vj  ;  •  «l  a+^jfKf-rtaHfts*- 


;rj-  - *r 
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For  the  camera  axis,  o  *  0.  Therefore,  its  azimuth  and  elevation  are 

-cob  +'c  Bin  (*g-*c) 

*“Ac  *  cos  ain#^  -  ain#^  cob  cos  ftg-Xc>  ' 

sinE  =  sin 4'  sin 4'  +  cos#'  cos#'  cos  ft  -X)  . 
c  sc  sc  sc 

For  hand  computations  it  is  convenient  to  write  Eq.  (100)  in  the  form 


^-cosE  cos  a' 

r  sin#'  cos  ft  -X  ) 

8  8  C 

sin  #^  sin(^g->-c) 

-cos#£ 

*  cos  12  > 

cos E  sinA 

X 

-sin  ft-X) 

8  C 

cos  ftfl-Xc> 

0 

sin  n  sin$ 

sinE 

cos#'  coa  ft  - X  ) 

cos  #'  sin  ft  -X  ) 

sin#' 

sin  12  cos  $ 

^  8  8  c 

8  sc 

S) 

i 

where 


(101a) 

<i°iy 


cos  12  *  cos#^  cos®  +  sin#^  sino  cos  (tj  +  fie) , 

sin  12  cos 4s  sin#'  cos o  -  cos#'  sino  cos  (rj  +  «e)  , 
c  c 

sin  12  sin  4  *  sin  o  sin  (ij  +  6e)  . 


It  can  then  be  readily  verified  that 

sinE  *  coso  sinEc  +  sino  cosEc  cos  (r,  +  6e  +  #) , 

tan  ft  -  A)  - - .  -#> - 7 - r 

c  coso  cos E^  -  sino  sinEc  cos(n  +  fie  -40  * 

where 


♦  =  tan 


‘h 


-cos  #'  sin  ft  -X  ) 
_S _ -J _ £_ 


sin  #'  cos  #'  -  cos  #'  sin  #'  cos  ft-  -  X 
sc  sc  sc 


J- 


(103) 

(104) 


(105) 


7.7  Horizon  to  Camera  Coordinates 

The  transformation  from  horizon  to  camera  coordinates  is  the  inverse  of  the 
transformation  given  by  Eq.  (100),  that  is. 
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/ -cosE  cosAV 

cosE  sinA  1.  (106) 


'V 


sinE 


Since  ordinarily  the  angular  rotation  of  die  film  plane  about  die  optic  axis  is  known, 
it  is  convenient  to  premultiply  both  sides  of  Eq.  (106)  by  R~*  (6e)  and  write 


{  sin  a  cos  (ij  +  6e)\ 
sin  a  sin(ij  +  6e) 


/  ®in^c  0  "C08*c\/C08<Xs'Xc>  -8to<VXc)  °\ 

sin  (X  -X  )  cos  (X  -X  )  0 

sc  sc 


Tc 

0  1 


\  COSO  /  \C08  4p 

0  sin  i'c/\  0  0 

/  sin  4’  0  cos4'\ 

s  S  | 

/-cosE  cos  A\ 

X 

0  10 

I  cosE  sin  A 

^  -cos  4'  0  sin  4'  ) 

8  S 

\  sinE  ) 

l 


(107) 


The  zenith  corresponds  to  a  direction  in  space  characterized  by  the  camera  co¬ 
ordinates  n  and  a  obtained  from  Eq.  (107)  by  setting  E  =  90°.  It  follows,  then, 

Z  Z 

that 


cos  a  *  sin  6'  sin  4'  +  cos  6'  cos  6'  cos  (X  -X  )  =  sinE_  , 

Z  CS  C  8  8  C  C 

-cos  $'B  sin(Xa-Xc) 

tan  +4c)  _  3iji6'  cos  4'  -  cos  41  sin  6'  cos(X  -X  )  * 


(108) 

(109) 


s 


8 


8  C 


In  Eq.  ( 108)  we  have  the  analytic  verification  of  the  obvious  relation  az  +  Ec  =  90  . 
For  the  general  case  we  can  write 


cos  a  -  sinE  sinE„  +  cos  E  cos  E  cos  (A  -  A  ) 
c  c  c 


(110) 


and 


tan  (q  +  fie-'l')  = 


cosE  sin(Ac  -  A) 


sinE  cosE„  -  cosE  sinE„  cos  (A  -  A) 
c  c  c 


(111) 


Comparison  of  Eqs.  (105)  and  (109)  shows  that  Eq.  (Ill)  can  be  written  in  the 
equivalent  form 


cos  E  sin  (Ac  -  A) 

tan(nz  -  v)  -  gjjjjj  COSE  -  cosE  sinE  cos  (A  -  A'  * 


(112) 


Explicit  expressions  such  as  Eqs.  (110)  and  (112)  are  useful  when  it  is 
desired  to  do  hand  calculations.  If  computerized  calculations  are  contemplated, 
it  is  more  efficient  and  straightforward  to  use  matrix  relations  such  as  Eq.  (106). 


This  is  because  we  can  have  computer  subroutines  for  the  general  rotations  Ry(6) 
and  R  (J),  and  can  generate  the  coordinate  transformations  by  successive  calls 


to  the  subroutines  with  appropriate  values  for  the  rotation  angles. 
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Appendix  A 

PmoMtiic  Cubic  Splinus 


In  the  ordinary  treatment  of  cubic  splines  (Ahlberg,  Nilson,  and  Walsh),  we 
consider  an  interval  a  £  x  <,  b  and  subdivide  it  by  a  mesh  of  points  corresponding 
to  the  locations  of  the  data  points: 


a  = 


X1<X2<‘ 


•<  xN+l  =  b* 


The  associated  set  of  ordinates  is 


yl*  y2’  y3*  •••'  yN+l  ’ 

We  denote  the  set  of  mesh  points  by  A  and  the  set  of  ordinates  by  Y.  We  then  seek 
a  function  S^(Y,  x)  that  is  continuous  together  with  its  first  and  second  derivatives 
in  the  interval  (a,b)  and  is  represented  in  each  subinterval  by  a  cubic  polynomial 
in  x.  If  we  write  the  Jth  polynomial  in  the  form 

y  -  y.  =  a^'x  -  x.)  +  b.(x  -  x^2  +  c.(x  -  x.)3  ,  i  *  1, 2,  3, ,  N  ,  (Al) 

the  spline  function  is  defined  when  all  the  coefficients  a.,  b^,  c^  are  known.  Since 
there  are  N  intervals,  the  total  number  of  unknown  coefficients  is  3N.  To  deter¬ 
mine  these  coefficients  we  have  3(N-1)  continuity  conditions  at  N-l  internal  mesh 
points.  A  further  condition  arises  from  the  Nth  polynomial,  which,  evaluated  at 

x  =  XN+1*  muSt  lead  to  y  =  yN+l' 
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A  unique  solution  requires  the  introduction  of  two  additional  conditions.  One 
of  the  simplest  ways  of  accomplishing  this  is  to  take 

sjP}  (a+)  *  S^(b-),  p  =  0, 1, 2  . 

The  spline  is  then  said  to  be  periodic,  of  period  (b-a).  For  nonperiodic  splines  the 
two  conditions  may  consist  in  specifying  the  end  slopes  or  requiring  that  the 
curvature  at  the  ends  be  zero.  It  can  be  shown  that  die  periodic  cubic  spline  with 
prescribed  ordinates  at  mesh  points  always  exists  and  is  unique.  The  same  is  true 
of  nonperiodic  splines  having  the  end  slopes  prescribed,  of  splines  having  prescribed 
second  derivatives  at  the  ends,  or  those  satisfying  the  conditions  S^(a)  =XS^(x2>; 
S^(xN+1>  =  p  S^(Xjj),  when  both  X  and  p  are  numbers  between  0  and  1. 

Spline  functions  composed  of  polynomials  of  type  (Al)  are  adequate  for  approxi¬ 
mating  single-valued  functions  of  x.  Multivalued  functions  can  be  handled  if  a 
parametric  representation  is  introduced.  The  most  suitable  parameter  is  die 
cumulative  chordal  distance  that  at  point  Pjkj*  Yj)  is 

j 

Sj  *  S  [fci  '  Xi-1)2  +  "  yi-l)2]l/2  *  3  *  2'3*4 . N+1  '  (A2) 

i=2 

The  problem  is  to  spline-fit  x  and  y  as  functions  of  the  parameter  s.  Note  that  the 
inclusion  of  multivalued  functions  requires  us  to  compute  two  sets  of  3N  polynomial 
coefficients. 

Let  us  denote  the  spline  functions  for  x  and  y  by  S^(X,s)  and  S^(Y,  s), 
respectively.  The  mesh  A  corresponds  to  the  set  of  s  values 

0  =  S1  <  s2  <  s3  <  '  *  •  <  SN+1  ’ 

The  cubic  polynomials  for  the  j.th  subinterval  are  then 


x  -  X.  =  afX\s-s.)  +  b^(s-s.)2  +  cfx\s-s.)3  , 

(A3) 

y  -  y.  =  af^fs-Sj)  +  bj^fe-s^2  +  cfy\s-Sj)3  . 

(A4) 

Because  of  the  relations 
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S-f  • 


dy  =  xjf  -jx  ' 
dx2  Ci)3 


where  the  dot  represents  differentiation  with  respect  to  s,  the  conditions  for 
matching  slopes  and  second  derivatives  for  the  function  y  =  y(x)  at  the  interval  mesh 
points  P^,  y.)  are  translated  into  corresponding  matches  of  the  slopes  and  second 
derivatives  for  each  of  the  spline  functions  S^(X,  s)  aud  S^(Y.  s). 

If  the  curvature  at  the  ends  is  to  vanish,  the  following  conditions  have  to  be 
satisfied: 


At  s  j,  x  =  y  =  0.  Consequently, 

blw  =  b Jrt  =  0  . 

At  sN+1,  x  ®  y  =  C.  Hence, 
bN)+3cN)(sN+l-8N)  =  0* 


b(y)+  3c(y)(g  -  s  )  =  o 

bN  N  'SN+1  SN'  0 


The  coefficients  for  the  x  and  y  splines  are  determined  by  repeated  application  of 
the  following  basic  procedure. 

Let  dt  =  si+1  -  st  and  let  be  either  (xi+J  -  x^/dj  or  (y.+1  -  y1>/di  .  At 
internal  points  we  have: 

f.  =  a.  +  d.b.  +  c.d.2  ,  i  =  1,2,. . . ,N  ,  (A5) 

ai+1  =  a.  +  2b.d.  +  3  c.d.2  ,  i  =  1, 2, . . . ,  N-l  ,  (A6) 

bi+1  =  bi  +  3c.di  ,  i  =  1,2,..., N-l  .  (A7) 

From  these  equations  it  follows  that 

dibi  =  3?i  -  ai+l  •  2ai  * 


so 
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This  tridiagonal  system  has  been  discussed  in  detail  in  the  literature  (Smirnov, 
1961).  The  recurrence  relations  that  are  obtained  are  well  adapted  to  digital 
computer  implementation. 

Figure  A1  shows  the  results  of  a  parametric  fit  to  the  synthetic  data  listed  in 
the  upper  right  corner  of  the  figure. 


Figure  Al.  Example  of  Curve-fitting  by  Means  of  Parametric  Cubic  Splines 
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Appendix  B 

Lagrange's  Expansion 

1-et  f(z)  and  g(z)  be  functions  of  z,  analytic  on  and  inside  a  contour  C  that  sur¬ 
rounds  the  point  z  *  a.  Let  €  be  such  that  the  following  inequality  is  satisfied  for 
all  points  on  the  perimeter  of  C: 

|€g(z)|<|z-n|  . 

Then 

u  =  a  +  eg(u) 

t  J  » 

has  precisely  one  root*  in  the  interior  of  C,  u  *  u(o);  and  f(u),  which  is  analytic 
on  and  inside  C,  can  be  written  in  the  form 

oo  n  ■  n-1  •  , 

f(u>  =  f(«)  +£  fr  [&)n  •  (BD 

,  1  ' 

This  is  the  Lagrange  expansion  (Bellman,  1964). 

*  This ’s  a  consequence  of  the  following  lemma:  Let  f(z>  and  <j>(z)  be  analytic 
in  the  interior  of  a  closed  curve  C,  continuous  on  the  curve  itself  and  such  that  on 
the  entire  curve  |<£(z)|<  jf(z)|  .  Then  the  two  equations  f(z)  =  0  and  f(z)  +  4>'z)  =  0 
have  the  same  number  of  roots  in  the  interior  cf  C, 
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Let  us  apply  this  result  to  Eq.  (86) 

2 

tant«(J)tan#,+^  Bint. 

with  e  changed  to  to  avoid  confusion  with  the  e  that  appears  in  Eq.  (Bl). 
We  introduce  the  following  identification: 

u  =  tant 

«“•(?«»  ■(?)— 7=  -(f) — Sj== 

cos  4'\1  +  tan  t  cos  4'^[ 1  +  u 

2 

c”€o 

f(u)  *  u,  f'(u)  *  1  . 

According  to  Eq.  (Bl)  we  can  then  write: 


-tsitn— -Hf 

1  \P  cos  4'  +  or  f 


.  e  /  ag 

*a  +  7 - 

\  p  cos  4 


■  fTs)*”  (»-♦')  (a ♦  «2)2’)+0<'4> ' 


that  is. 


tant-(5)tan*'+«*(&Xr)  p=^=  -  + 

l1  ■  «o  sin  *7 


Substitution  of  ~  yl-e^  leads  to  Eq.  (87)  . 


Appendix  C 


Useful  Constants 


Ti«e  Coostaats 

Period  of  rotation  of  earth:  23h56m04.  Is  mean  Bolar  time 

One  uniform  sidereal  day:  23h56m04. 0918  mean  solar  time 

One  mean  solar  day:  24h03m56.  555s  uniform  sidereal  time 

Tropical  Year:  366.2422  mean  sidereal  days 
365. 2424  mean  solar  days 

Sidereal  Year:  366. 2564  mean  sidereal  days 
365, 2504  mean  solar  days 


Astrodynaaic  Constants 


Rotational  rate  (w)  *  0.729211514  (6)  X  10"4  rad/sec 
Mass  of  earth  (M)  *  5. 975(1±0. 0007)  x  1027  gm 
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Constants  (or  Geographic  Ellipsoids 


Author,  Date 

ae(m) 

1/f 

Countries  using  ellipsoid 

Clarke.  1866 

6,378,206 

294.98 

United  States,  Canada, 
Mexico 

Clarke,  1880 

6,378,249 

293.47 

France,  South  Africa 

Everest,  1830 

6,377.253 

300.  80 

India 

Plessis 

6,376,523 

308.  64 

France  (mapping) 

Bessel,  1841 

6,  377, 397 

299.  15 

Germany,  Austria,  Dutch 
East  Indies 

Kraijenhoff 

6.376, 950 

309.  65 

Holland 

Danish  survey 

6,377,019 

300 

Denmark 

Constants  (or  International  Ellipsoid 

Semimajor  axis  (equatorial  radius)  *  6,378,388  m 

Flattening  =  1/297  =  0,003  367  003  4 

Polar  radius  =  6,356,911.946  m 

Square  of  eccentricity  =  0.006  722  670 

Length  of  quadrant  of  the  equator  =  10,019, 148.  4  m 

Length  of  quadrant  of  the  meridian  *  10, 002, 288.  3  m 

Area  of  the  ellipsoid  *  510, 100,934  sq  km 

Volume  of  the  ellipsoid  =  1,083,319,780,000  cu  km 

Radius  of  sphere  having  same  area  as  ellipsoid  *  6,371,227.7  m 

Radius  of  sphere  having  same  volume  as  ellipsoid  =  6, 371,221. 3  m 

2 1 

Mass  of  the  ellipsoid  =  5.  998  x  10  metric  tons 
Mean  density  =  5.  527 


